Solving the Power Flow Problem in Bipolar DC Asymmetric Distribution Networks Using Broyden’s Method

This research addresses the power flow analysis in bipolar asymmetric direct current (DC) networks by applying Broyden’s numerical method. This general successive approximations method allows for a simple Newton-based recursive formula to reach the roots of multiple nonlinear equations. The main advantage of Broyden’s approach is its simple but efficient structure which can be applied to real complex nonlinear equations.The power flow problem in bipolar DC networks is still challenging, as multiple operating options must be considered, e.g., the possibility of having a solidly grounded or floating neutral wire. The main goal of this research is to contribute with a generalization of Broyden’s method for the power flow solution in bipolar DC networks, with the main advantage that, under well-defined conditions, this is a numerical method equivalent to the matricial backward/forward power flow, which is equivalent to the successive approximations power flow method. Numerical results in the 21-, 33-, and 85-bus grids while considering two connections for the neutral wire (i.e., solidly grounded at any node or floating) show the effectiveness of Broyden’s method in the power flow solution for bipolar asymmetric DC networks. All numerical simulations were carried out in the MATLAB programming environment.


General Context
Bipolar DC networks are emerging technologies that aid in providing electricity in medium-and low-voltage networks by using three wires associated with the positive (p), neutral (o), and negative (n) poles [1]. These networks can be considered to be the DC equivalent of the conventional three-phase AC distribution networks given that multiple constant power loads can be connected between their poles. These may be monopolar (connected between the positive or negative pole and the neutral wire) or bipolar (connected between the positive and negative poles) [2,3]. Figure 1 illustrates the typical configuration of a bipolar DC network with multiple constant power terminals [4].
The main characteristics shown in Figure 1 are: i. The operation of a bipolar DC network can be unbalanced due to the presence of constant monopolar and bipolar power terminals, which implies that no reductions (monopolar equivalents) must be applied [5]. ii. Concerning purely monopolar configurations, bipolar DC grids imply an increase of about 34% in construction investments. However, more than two times the power can be transferred by these grids, including the possible connection of exceptional loads, i.e., bipolar compositions [2].
iii. Depending on the operating practices of the distribution company, the neutral wire can be operated in two main ways: floating in all nodes except the substation bus or solidly grounded at any network node [6].

Motivation
Considering the main characteristics of bipolar DC networks, including their unbalanced operation and the possible connections of the neutral wire, power flow studies for multipolar DC networks with constant power terminals can be modeled as a set of nonlinear equations that must be solved through the application of numerical methods, as their analytical solution is impossible [7]. This research is motivated by applying Broyden's (also known as the secant) method for finding roots in simultaneous sets of nonlinear equations [8]. The selection of this method is based on its easy formulation and implementation, as well as on its linear convergence, with the main advantage that infinite possible initial matrices A t (with t = 0) may be selected for its numerical validation.
The authors of [11] analyzed a power flow method under different operation modes and control strategies in bipolar high-voltage DC (HVDC) grids integrated with voltage source converters (VSCs). The power flow method was based on layers and the admittance matrix of the bipolar grid.
The authors of study [7] proposed a successive approximations method to solve power flow in bipolar DC distribution networks. This method allows considering a grounded or non-grounded neutral wire.
The authors of study [3] demonstrated the exact conditions and necessities for the convergence and uniqueness of the power analysis solution in bipolar DC grids. This study also considered the over-current protection scheme of power electronic converters.
The authors of work [6] tackled the power flow problem in bipolar DC grids with multiple constant power loads by employing a triangular-based power flow formulation. This formulation is suitable for radial setups and considers whether the neutral wire is grounded or not at all system nodes.
The authors of study [7] used a successive approximations technique based on the classical backward/forward method to solve the power flow analysis in bipolar DC grids, including the possibility of grounding or not grounding the neutral wire. Furthermore, the convergence of the proposed approximation was demonstrated using the Banach fixed-point theorem.
The authors of [12] presented a formulation to solve the power flow problem in bipolar DC networks, which used the Taylor expansion in the bus injection equations, generating a hyperbolic formulation that can be implemented in both meshed and radial networks.
Previous studies have shown promising results, but some of them struggle to achieve convergence in their performance. Additionally, some of these studies lack a comprehensive framework to address these limitations. This paper presents a solution to the power flow problem in bipolar DC asymmetric distribution networks using Broyden's method. This method is a quasi-Newton method that numerically solves systems of nonlinear equations through iterative processes. In addition, it offers several advantages over other approaches [13,14]. For example, this method provides greater versatility, making it applicable to a broader range of problems without requiring additional modifications. Additionally, Broyden's method has demonstrated good and rapid convergence. Moreover, it is robust, as it can be implemented in cases where other techniques, such as Newton's method, may encounter difficulties due to ill-conditioned Jacobian matrices.

Contributions and Scope
Considering the state of the art in numerical methods applied to the solution of the power flow problem in bipolar DC networks with asymmetric loading, the main contributions of this research are the following: i.
The application of Broyden's method to the power flow problem in bipolar DC grids with the aim to solve the set of Equations ∇ f (x) = 0 by using different A t matrices. ii. The evaluation of the convergence properties of Broyden's method using random initialization matrices A t with t = 0, as well as an analysis of its equivalence with the successive approximations method when A 0 is selected as part of the conductance matrix of the system. iii. A demonstration of the equivalence between the successive approximations power flow (SAPF) and Broyden's method, which proves that the former is a particular case of the latter in the studied problem.
Regarding the scope of our contribution, the following facts are considered during the numerical validations carried out in the 21-, 33-, and 85-bus grids: (i) two possible connections for the neutral wire are tested: the first case corresponds to the floating connection, i.e., the neutral wire is only solidly grounded at the substation terminals, and the second case assumes that all the nodes are solidly grounded; (ii) the load consumption values are assumed to be perfectly known, i.e., no uncertainties are considered in the values of the constant power terminals; and (iii) the voltage output at the substation bus for all the poles is fixed as an input parameter.

Document Structure
The remainder of this document is structured as follows. Section 2 presents the general formulation of the power flow problem for bipolar DC networks with asymmetric loading. Section 3 describes the general implementation of Broyden's method for a general set of nonlinear equations. Section 4 presents the application of Broyden's method to the power flow problem. Section 5 shows the main characteristics of the test feeders, which correspond to the 21-, 33-, and 85-bus grids. Section 6 presents the numerical validations of the proposed power flow problem, including the voltage profile performance and the power losses for each test feeder, as well as some numerical comparisons against the literature reports. Finally, Section 7 lists the main concluding remarks derived from this work and some possible future studies.

Power Flow Problem in Bipolar DC Networks
The power flow solution for bipolar DC networks with asymmetric loading is a challenge due to the presence of multiple nonlinear equality constraints associated with each pole as a function of the connection of multiple constant power loads (monopolar or bipolar) [15]. This set of constraints requires implementing numerical methods in order to obtain a solution with an acceptable convergence error [1]. The power flow problem in a bipolar unbalanced network can be represented with seven equations (three linear and four nonlinear) per node and pole. These equations are presented below.
where I p g,k , I o g,k , and I n g,k represent the current injection in the slack source per pole (i.e., positive, neutral, and negative -pon-); I p d,k , I o d,k , and I n d,k are the current demanded by the monopolar constant power loads, i.e., P p d,k and P n d,k (note that the monopolar loads refer to those connected between the positive or negative pole and the neutral pole); I p−n d,k corresponds to the current absorbed by the constant power load P p−n d,k connected between the positive and negative poles, i.e., a bipolar consumption; I ground d,k means the current drained from the neutral wire to the physical ground at each node; V r k represents the voltage magnitude at node k in pole r; G pr jk , G or jk , and G nr jk correspond to the conductance values (obtained from the conductance matrix) that relate nodes j and k and poles pon with pole r, respectively. Note that P and N denote the sets that contain all the poles and nodes of the bipolar DC grid under analysis.
Remark 1. The set of nonlinear Equations (8)-(10) can be compacted using a matricial representation with the help of the conductance matrix and its components G gg , G gd , G dg , and G dd , as defined in (11) [7].
where, considering that the cardinalities of the sets P and N are 3 and n, I g ∈ R 3×1 and V g ∈ R 3×1 are the vectors that contain the current injection in the slack source (unknown variables) and the voltage output at the terminals of the slack source, i.e., perfectly known voltages.
is the vector containing the current demanded at each node, which is ordered by node and pole (unknown variables); and V d ∈ R 3(n−1)×1 is the vector of voltage profiles at each demand node, also ordered by node and pole (unknown variables).
From the two equations in (11), it can be noted that: i. The first row of (11) shows that the solution of the current injections at the substation will be only known when all the demanded voltage profiles are determined. In addition, this is a linear equation that does not require any iterative process for its solution. ii. The second row of (11) is a nonlinear set of equations, as the demanded currents I d are hyperbolic functions of the demanded voltages, which implies that its solution requires efficient numerical methods [1].
To determine the set of demanded currents by comparing the second row of (11) against Equations (8)-(10), note that where I B d and I M d are vectors associated with the currents demanded by the bipolar and monopolar loads, respectively, which are calculated per node, as follows: where

Remark 2.
The solution to the power flow problem in bipolar DC grids with asymmetric loading corresponds to the solution of the following recursive formula: using a numerical method based on the dependence of the vectors I B d and I M d on the demanded voltages V d .

The Secant or Broyden's Method
The secant or Broyden's method is a numerical approach developed for solving sets of nonlinear equations with ∇ f (x) = 0. This is performed without resorting to the implementation of Newton-based approaches [8,16]. The general successive approximations formula for the secant method is defined in (6).
where t is the iteration counter and A t corresponds to the Jacobian matrix, which is, however, approximated via secant hyperplanes [17]. Consider the following linear relation: which is simplified as (17) by Now, the following two hyperplanes are defined: The main idea of the secant method is that, at the solution point, hyperplanes (18) and (19) must be equal, i.e., Note that, in (20), if x = x t is defined while considering the definition in (17) and making some algebraic manipulations, the following result is reached: By simplifying (21), the general updating formula for A t can be determined as presented in (22).
Note that, for the general iteration t + 1, the way to update A t+1 is obtained from (22), as defined in (23).
The general algorithm for solving the set of equations ∇ f (x) = 0 is presented in Algorithm 1 [18].

Algorithm 1: General implementation of Broyden's method for solving sets of nonlinear equations
Data: The set of nonlinear equations under analysis, i.e., ∇ f (x) = 0. Result: The solution vector x * . 1 Make t = 0 and define an initial estimation of A 0 ; 2 Define the maximum number of iterations t max and the convergence error ε;

Application of Broyden's Method to the Power Flow Problem
This section presents the general application of Broyden's method to the problem regarding the power flow of bipolar DC networks with asymmetric loading, which is based on Algorithm 1. As presented in Section 2, this work deals with a complex nonlinear problem, where the main idea is to efficiently obtain its roots via numerical methods [1]. For the application of Broyden's method, as described in Section 3, using the compact representation in Equation (14) is recommended. Note that this Equation is necessary as the equivalent gradient function ∇ f V t d is defined from it, as follows [18]: where the components of the demanded currents in I t,B d and I t,M d are defined from (12) to (13).
Considering the structure of Algorithm 1, Algorithm 2 presents the application of Broyden's method to solve the power flow problem in asymmetric bipolar DC distribution grids.

Remark 3.
The main characteristic of Algorithm 2 is the selection of the initial value of matrix A 0 , as it has multiple alternatives. Depending on this selection, the number of iterations required to reach the desired convergence varies. Obtain the components G dg and G dd from G bus ; . . , n, k = s; 7 Define the maximum number of iterations t max ; 8 Define the convergence error ε; 9 for t = 0 : t max do 10 Construct the current vectors in (14) by using (12) and (3) ∀k = 1, 2, . . . , n, k = s;

Test Feeders
To validate the effectiveness of Broyden's method in addressing the power flow problem in bipolar asymmetric distribution grids, three test feeders composed of 21, 33, and 85 nodes with unbalanced structures were considered.

Bipolar DC 21-Bus Network
The 21-bus grid corresponds to a radial distribution network composed of 21 buses and 20 lines, as depicted in Figure 2 [6]. The main characteristics of this feeder are the following: (i) the substation bus is located at bus 1, and the voltage outputs in the positive and negative poles are ±1 kV with respect to the neutral point, which is solidly grounded; and (ii) the total monopolar consumption is 554 kW in the positive pole and 445 kW in the negative pole, while the bipolar power consumption is 405 kW.
The information regarding branches and monopolar and bipolar power consumptions is presented in Table 1.

Bipolar DC 33-Bus System
This bipolar DC network is a modification of the original single-phase distribution network originally proposed by the authors of [19]. The schematic single-line representation of this system is presented in Figure 3. The main characteristics of this distribution network are the following: (i) the substation bus is located at bus 1, and the voltage outputs in the positive and negative poles are ±12.66 kV with respect to the neutral point, which is solidly grounded; and (ii) the total monopolar consumption is 2615 kW in the positive pole and 2185 kW in the negative pole, while the bipolar power consumption is 2350 kW.   Table 2 [20].

Bipolar DC 85-Bus Network
The bipolar DC 85-network is a radial distribution network composed of 85 nodes and 84 branches, as depicted in Figure 4. The main characteristics of this distribution network are the following: (i) the substation bus is located at bus 1, and the voltage outputs in the positive and negative poles are ±11 kV with respect to the neutral point, which is solidly grounded; and (ii) the total monopolar consumption is 1745.48 kW in the positive pole and 2682.19 kW in the negative pole, while the bipolar power consumption is 2258.58 kW.
The information regarding branches and monopolar and bipolar power consumptions is presented in Table 3 [21].

Numerical Validations
To validate Broyden's method, all the bipolar DC networks were simulated in the MATLAB programming environment via the researchers' own scripts. For this implementation, version 2021b was used on a PC with an AMD Ryzen 7 3700 2.3-GHz processor and 16.0 GB RAM, running on a 64-bit version of Microsoft Windows 10 Single language.

Number of Iterations and Convergence Behavior
In this section, the A 0 matrix is generated using the conductance matrix G dd with the following rule: (1) is a random number generated using a uniform distribution, and the parameters α and β are set as 0.5 and 2.0, respectively. In addition, to test the variations in the number of iterations required to reach the power flow solution while considering a convergence error of about ε = 1 × 10 −10 , for these simulations, it is assumed that the neutral wire is floating in all nodes except in the substation bus. The number of iterations required for solving the power flow problem with Broyden's method in each one of the test feeders is depicted in Figure 5. The main characteristics of the behavior shown in Figure 5 allow stating that: i.
The number of iterations after 100,000 repetitions exhibits a Gaussian distribution, with a mean value of about 11 iterations for all the test feeders. This is due to the fact that the Gaussian distribution carried most of the γ-factor around its center, i.e., α+β ii. Depending on the γ-factor, for all the test feeders, a minimum number of iterations of about seven was observed. The average values of the γ-factor were 0.9750, 0.9659, and 0.9783, for the bipolar DC 21-, 33-, and 85-bus networks, respectively. These values imply that in order to reach a better numerical performance with Broyden's method, the γ-factor must be tuned for each test feeder. iii. The maximum number of iterations for the 21-bus grid was about 14, whereas for the 33-and 85-bus grids it was 13. For these results, the γ-factor is when this parameter is located near the α of β values, i.e., near the extreme values used in simulations.
To determine the convergence properties of Broyden's method in all the test feeders under analysis, a logarithmic error graph is presented in Figure 6. Equation (25) presents the function that defines the error's evolution.
( The convergence behavior exhibited by Broyden's method in Figure 6 for all the tested feeders when the γ-factor was selected to ensure the minimum or the maximum number of iterations within the range assigned for α and β shows that, in general, this algorithm has a linear tendency to solve the power flow problem, with some oscillations. However, these oscillations are caused by the updating of the A t matrix, which is carried out by using the information of the previous iteration. Still, this does not compromise the stability of this method.

Remark 4.
The γ-factor is a key parameter in the implementation of Broyden's method. However, as observed for all the test feeders, the recommended value to ensure the minimum number of iterations must be contained in the 0.90-1.10 interval.

Comparative Analysis
To demonstrate the effectiveness of the proposed secant method to deal with the power flow problem in bipolar DC asymmetric distribution networks, the power flow was solved while considering the two connection modes of the neutral wire (i.e., solidly grounded or floating) in the tested feeders. A comparison was made with the hyperbolic approximations power flow (HAPF) approach reported in [6], the successive approximations power flow (SAPF) method proposed by [7], and the triangular-based power flow (TBPF) approach reported in [12]. Table 4 lists the results for all the test feeders, including the two connections applicable to the neutral wire. For the sake of simplicity, the gamma factor was selected to ensure minimal iterations in Broyden's power flow method. The numerical results in Table 4 show that: i. In all simulation cases, Broyden's method reaches the same numerical solution regarding power losses in the tested feeders for both neutral wire operating conditions. Nevertheless, when the neutral wire is floating, this approach takes 7 iterations; otherwise, 11 iterations are needed. These values were independent of the number of nodes of the test feeder under analysis. ii. Contrary to the behavior of the comparison methods, the number of iterations increases when the neutral wire is considered to be solidly grounded. However, this behavior can be attributed to the fact that Broyden's recursive formula is general and does not depend on the nonlinear set of equations under analysis, which implies that, under some particular conditions (i.e., the connection of the neutral wire), its evolution differs from that of specialized methods for power flow studies. On the other hand, the increase in the number of iterations is not directly related to that of the required processing times. iii. The main characteristic of the simulations in the three test feeders is that (as expected) when the neutral wire is solidly grounded, the power losses are lower than in the floating operation scenario. These differences are 4.1536, 10.0629, and 37.2778 kW, for the 21-, 33-, and 85-bus networks, respectively.
On the other hand, the behavior of the voltage profiles in the tested feeders shows voltage regulations of 11.1740%, 9.4264%, and 7.9629% for the neutral floating case in the 21-, 33-, and 85-bus networks, as well as 10.9898%, 9.3176%, and 8.1086% when the neutral wire is solidly grounded. As expected, voltage regulation improves when the neutral wire is solidly grounded, which is a direct consequence of the reduction in the number of power losses evidenced in Table 4. i.
The number of iterations required by the proposed numerical method to solve the studied problem highly depends on the γ factor selection. This parameter must be tuned for each test feeder and neutral wire connection, with its recommended values being between 0.90 and 1.10. Note that the γ-coefficient can be considered to be the equivalent value of the α-coefficient in the classical Gauss-Seidel method implemented for AC power systems. ii. The convergence behavior, i.e., the evolution of e(t), exhibits a linear behavior as a function of the selected γ-factor, with slight oscillations attributed to the updating rule applicable to Broyden's method, which adds the effect of the difference between two consecutive nonlinear function evaluations to the calculation of A t . In addition, after 100 thousand consecutive evaluations of Broyden's method for different values of the γfactor, it was observed that the highest probability in the expected number of iterations was 11, regardless of the test feeder, with a minimum value 7 for all test systems, as well as a maximum of 14 for the 21-bus grid and 13 for the remaining grids. iii. Numerical comparisons with the literature reports showed that the processing times in the solution of the power flow increase with the number of nodes of the network, as expected. However, regarding the number of iterations, the opposite behavior was evidenced in the methods used for comparison (SAPF, TBPF, and HAPF): in the case of the floating neutral wire, the number of iterations was lower for our method, and, when the neutral wire was solidly grounded, the number of iterations increased. This, however, was inversely related to the average processing times. iv. A demonstration of the equivalence between the SAPF approach and Broyden's method under Assumptions 1 and 2 confirmed the generality of the proposed technique in dealing with power flow problems in bipolar DC networks with asymmetric loads, with the main advantage that the number of iterations or processing times can be prioritized as a function of the γ-factor and the selection and adapting of the A t+1 matrix.
The main advantage of Broyden's method for addressing the power flow problem in bipolar DC networks with asymmetric loading is that it can apply to radial and meshed networks as an equivalent to the SAPF and the HAPF approaches given that its formulation is based on a conductance matrix that is constructed without considering the tree configuration of the network under analysis. However, numerical tests were only conducted in three radial distribution networks since these configurations correspond to the worst possible case regarding power losses and voltage regulation indicators.
In regard to future works, the following studies can be conducted: (i) extending the proposed power flow method to three-phase asymmetric distribution networks; (ii) including, in Broyden's formulation, the possibility of having multiple voltage-controlled nodes and distributed energy resources; and (iii) demonstrating the equivalence between the SAPF method and the proposed approach. Funding: This research received support from the Ibero-American Science and Technology Development Program (CYTED) through thematic network 723RT0150, Red para la integración a gran escala de energías renovables en sistemas eléctricos (RIBIERSE-CYTED).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing does not apply to this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
Indices p, o, n Superscripts associated with poles 0 Superscripts associated with the initial value t Superscripts associated with the iteration counter j, k Subscripts associated with nodes g Subscripts associated with the generators d Subscripts associated with the demands Parameters Vector that contains the current in monopolar loads (A) A t Jacobian matrix in iteration t